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Abstract 

Lowest-order cumulants provide important information on the shape of the emission 
source in femtoscopy. For the simple case of noninteracting identical particles, we show how 
the fourth-order source cumulant can be determined from measured cumulants in momen- 
tum space. The textbook Gram-Charlier series is found to be highly inaccurate, while the 
related Edgeworth series provides increasingly accurate estimates. Ordering of terms com- 
patible with the Central Limit Theorem appears to play a crucial role even for nongaussian 
distributions. 
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1 Introduction 

The large experimental statistics which are now available permit femtoscopic correlations of 
identical particles (see e.g. Ref.pQ) as a function of the full three-dimensional momentum dif- 
ference q = pi — P2 and often also of the average pair momentum K = (pi + P2)/2. Increasing 
attention has therefore been paid to the detailed description in these higher-dimensional spaces 
of the second-order correlation function 



p(q,K) 

P rcf (q,K)' 



l + fl(q,K) = C(q,K) = ™ ^ , (1) 



with p the density of like-sign pairs in sibling events and p rei the reference pair density usually 
determined by a combination of event mixing and Monte Carlo simulation. After removing 
irrelevant correlations, [2] the correlation function can yield information on the spacetime sta- 
tistical properties of particle emission as embodied in the source function S(x, K). The source 
function is obtained from the two-particle emission function, the density of particle emission 
points with near-equal momenta, by projection from the two four-coordinates onto the relative 
three-coordinate x (as measured in the pair rest frame of emission point densities) [3j SI [5j [6] . 
The momentum- and coordinate-space correlations of noninteracting identical particles in their 
centre-of-mass system are related by a Fourier transform 



R(q,K) = X J d 3 iS(x,K) e* 
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where A is the correlation strength parameteiQ. Suppressing the K-dependence in a given global 
system, the two-particle correlation in both spaces can in this case be written as a normalised 
probability density function (pdf) /(q) in g-space and a pdf g(x) in x-space, related by 



/(q) 



fl(q) 




(3) 




fd?qR(q) 
S(x) 



5(x) 



(4) 



/AS(x) 



(2tt) 3 S(0)/(0)- 



A gaussian /(q) immediately yields a gaussian g(x) in any dimension. Experimental data, 
however, is often nongaussian, sometimes strongly so. This raises two problems: first, to sys- 
tematically describe the nongaussian shape of R (or /) in momentum space, and second, to 
determine parameters of S (or g) in coordinate space, given only the kernel transform and 
measurements in g-space. 

Approaches towards systematic description of nongaussian shapes in q-space can be found 
in e.g. Refs. [II El E]) while the source function S is reconstructed by means of higher-order 
coefficients in x-space using imaging techniques [TOj, [UJ 02] and cartesian harmonics [21 [13] . 

In this paper, we wish to address the second problem of a systematic description of g(x) 
in terms of given measurements in g-space, based on the fundamental statistical properties of 
cumulants; the corresponding approach treating the first problem of measurements in g-space 
has been treated in part in the literature p3] [15] and will be more fully elaborated elsewhere. 

2 Cumulants in dual spaces 

While fully three-dimensional formulations have been in part set out in e.g. Ref. [15] . we shall here 
work in one dimension using so that the above expressions become g(x) = j dq e~ tqx f(q) /2nf(0) 
and so on, our purpose being first to test and improve the convergence properties of series 
expansions in a simpler environment. 

Given a measured normalised correlation function f(q), its g-moments //j- 9 ' = f dqf(q)q r 
and g-cumulants of lowest orders r = 1,2,3,... provide fundamental information on its 
properties: the ordinary mean fi^ = = J dq f(q) q is a measure of the location of the 
peak of f(q), while the variance = ^ ~ (Mi ) 2 measures the dispersion and a = (k^ 9 ') 1 / 2 
the width or scale of the pdf, the skewness = K^/a 3 measures its asymmetry and the 
kurtosis "f^ = /<t 4 is a first description of the pdf tail's decay rate. Higher-order "generalised 
kurtoses" 7r 9 ' = 8 /<r r would provide successively more detail. Kurtoses 7^ can also be 
generally viewed as cumulants of the pdf /(</) of the standardised variable q' = (q — /ij )/<r. 

Equivalent relations hold in coordinate space between x- moments, x-cumulants and g(x), 



e.g. fj,r x) = J dx g(x) x r , = /-4 — (a*i ) 2 and so on. 

g-moments are derivatives of the generating function $(a?) = 27r/(0)g(— x) = j dqe iqx f(q), 



For example, A = 1 and ( — 1/2) for noninteracting identical spin-0 and unpolarized spin- 1/2 particles emitted 
at independent spacetime points respectively. The A-parameter may be modified by particle impurity or a contri- 
bution of the particles emitted in some multiparticle quantum states analogous to the coherent states in quantum 
optics. Note that A is eliminated by the respective normalisation in Eqs. ©-((3J). 




(5) 
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writing D r x = {d/dx) r for short, while the related derivation of (/-cumulants from 

K^ = (-iYD r x ln^ X )\ x=0 (6) 

fixes relations between moments and cumulants to all orders. For identical particles, both 
C(q) = C(— q) and g(x) are symmetric, so that moments and cumulants of odd order vanish 
and the even-order relations in both (/-space and x-space become 

K2 = HQ, (7) 
K 4 = fj, A - 3f4, (8) 

^6 = — 15/^4^2 + 30/^2 etc. (9) 

Cumulants form a natural basis for near-gaussian expansions since a gaussian pdf is fully deter- 
mined once K\ and K2 are known: all its K r >3 are identically zero. They also have important 
properties such as invariance under translation and a null result for uncorrelated variables. 

While for purely gaussian sources, the second-order cumulants are related by k£ = l/^^ 
and all higher-order cumulants are identically zero, neither of these statements is true in general. 
We will therefore consider both the modification of « 2 resulting from nonzero 7r 9> as well as 
the x-kurtosis = «4 /( K 2 ) 2, Since x-moments are found from the generating function 

*(?) = /(?)//(0) (10) 

through 

^ = (-iYD^(q)\ q=0 , (11) 
we can through Eqs. (J7])-(|9]) obtain x-cumulants as combinations of measured (/-moments. 



3 Gram- Char lier expansions 

3.1 Expressing 7^' in terms of 7^' 

While experimental measurement of derivatives of &(q) is of course impossible, the above ex- 
pressions can nevertheless be evaluated since Gram-Charlier and Edgeworth series expansions 
also probe regions of nonzero q. Both expansions start with choosing a reference pdf fo(q) 
which, given the close relation between cumulants and gaussians, is almost invariably chosen by 
textbooks [161 E] to be a gaussian 

h (q) = e —^ i.e. fo(q') = — / =- (12) 



with the free parameter a 2 fixed to the experimentally measured ■ The resulting "Gauss 
Gram-Charlier" (GGC) series, also known as the "Gram-Charlier Type A" series, and the cor- 
responding Gauss Edgeworth (GEW) series are closely related, being mere re-orderings of one 
another, and are therefore commonly considered to be one and the same. As we will show, 
however, the GEW far outperforms the GGC series at any order of the partial sums. 
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As shown elsewhere. |18j the GGC series results from expanding the generating function for 
the nongaussian f(q') in powers of x' 



e x ' 2/2 exp 



j=3 



-x' 2 /2 (ix') m 



m=0 



(13) 



where each c m (£) is a polynomial in the set of g-kurtoses £ = {( r = 7r <?) / r '}?=4- Taking 
the inverse Fourier transform of &(x') term by term, one obtains an expansion in terms of 
Chebychev-Hermite polynomials H r (q'), 



/(«') = kW) 



i=2 



Qy(C) 
(2j)! 



flrCfl') = fo\<l') (-D q >) r fo(q') , 
with lowest-order terms (writing H r {q') = H r for short) 



(14) 
(15) 



fW) = foW) 



Using {-D ql yUq')H 2j (q') 
function is, for even r, 



1 + C 4 #4 + ( 6 H 6 + (Cs + ^(i)H 8 + (Cio + Ce Q)H 
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(16) 



fo(q')H2j+ r (q'), the r-th derivative of the x-moment generating 



= e - 9 ' 2 /2 ^ r( ^ } + C4jff4+r((? ') + C 6 # 6+r ( g ') + ...], 

from which the x-cumulants follow as ratios of generating functions at q' 
generalised g-kurtoses 7r = /a r and using i?2r(0) = (— l) r (2r — 1)!! 



in terms of 



(-i)2$(2)(g/) 



.(«) 



while the x-kurtosis in fourth order is 

(x) _ $(4) $ (0) _ 3$ (2) $ (2) 
74 ~~ $(2) $(2) 

with a similar expression for k£ . 

.,(*) 



1 



1 + §74 - ^76 + 1^(78 + 35 7 |) + . . . 
1 + 174 - ^76 + sk(78 + 35 7 |) + . . . 



74 



he + 578 + x74 



1 + 174 



2^76 + ^78 + §7 4 2 + • • 



(17) 



(18) 



Note firstly that depends only on 7^ 9) but not directly on a 



v 2 ' 



this is true also for 



higher-order . Secondly, the above relations reduce to the gaussian relation = 1/ K 2^ an d 
7 4 a:) = if and when the measured correlation function is gaussian since as mentioned all 7^> 3 

are then identically zero. In general, however, the "radius" [2k 2 J 1 / 2 of the source distribution 
is a function also of higher-order q-cumulants, with both increasing orders 7 r 9> and increasing 
powers of lower-order entering the expansions. 

Given the symmetry between x and q, the corresponding expansions for t6f and jr q) in terms 
of and 7^ would have the same form as the above, apart from some changes in sign. Any 
measured k£' is therefore itself the result of contributions from higher-order cumulants of g{x) 
or, in physics terms, the nongaussian shape of the emission region. 
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3.2 Truncation and the GGC disaster 



Statistical errors on g-cumulants rise with increasing order so that only those lower-order ones 
accessible to available experimental statistics can be included. Series expansions such as (114p 
are known to be asymptotic, so that the question arises: how accurately can a series truncated 
at some maximum order 7 rmax and/or a maximum power 7^ max estimate the jr ■ 

To quantify this issue, we make use of the Normal Inverse Gaussian (NIG) probability density 
|19j as a solvable toy model for f(q') which yields exact expressions for both coordinate- and 
momentum-space cumulants. While the NIG has four parameters a, /3, \i and 5, in the present 
symmetric case /3 = (J, = 0, so that we need only the two-parameter Symmetric Normal Inverse 
Gaussian (SNIG), 

f(q\oi,d) = , (19) 

ttVS 2 + q 2 

where K\ is the modified Bessel function. The SNIG reverts to a gaussian in the limit a — > oo and 
has g-moment generating function <!>(x | a, 5) = exp[5a — #V a 2 + x 2 ] . Experimentally measured 
and 74 ?) would fix the parameters: writing a = [k^] 1 / 2 and g = 7 4 9) for short, a = [3/ga 2 ] 1 / 2 
and S = [3a 2 / g} 1 / 2 , so that higher-order cumulants and kurtoses can be expressed in terms of 
measured quantities a and g as 



k$mg = (r-iy.\(r-3)\\a r [lg}2-\ (20) 
7$ MG = (r-l)!!(r-3)!! [U]^ 1 . (21) 



L3 

Using the SNIG pdf as x-moment generating function in the form ([10 



we obtain exact expressions for x-cumulants via (|lip and the moment-cumulant relations. Omit- 
ting the argument of the Bessel functions, which is ad = 3/g in every case, these "exact" 
a>cumulants are 



.to 



1 K 2 

' "2,SNIG " -^i'Kl* 

M J_ 3K Z K X -3K 2 

K 4,SNIG — (V>2 ' TS2 ' V ZO J 



,(8)2 K 2 
"2 1 

3K 3 Kij- 3K\ 

l 2 



7g NIG = ""°" K 2 — ■ ( 24 ) 



With these exact x-cumulants as reference, we test the accuracy of various truncations of 
Eqs. (|17[) - (|18p as a function of the Gram-Charlier order m = 2j of Eq. (|14p . 

The results are disastrous. In Fig. 1, we show respectively the percentage deviation of GGC 
expansions (|17|) and (|18p . truncated at mth order, from the exact answers (|22p and (|24p . in the 
form IOO^^/k^snig ~~ 1) anc ^ 100 (74 ^/ 74 snig ~~ -0- At 74 = °^ course ) an series reduce to a 
gaussian and all approximations become exact. Even small values of 74'' lead to large deviations, 
however, and the size of the deviations increases with order m. GGC series fail completely to 
approximate the exact x-cumulants. 
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Figure 1: Gauss Gram-Charlier series fail badly: Percentage deviations of Gram-Charlier ap- 
proximations of «2 found by Eq. (fTTj) vs. ([22]) and of ' found by (fl~8j) vs. (f24"|) . as a function 
of the measured g-kurtosis 7 4 9) for various Gram-Charlier orders m. 



4 Edgeworth series 

4.1 Derivation and properties 

In his 1946 treatise on statistics, Cramer [20J derived the Gauss-Edgeworth (GEW) series by 
considering the random variable (q' in our case) to be a convolution of n identically-distributed 
independent (iid) random variables qi each with pdf fi(qi), a corresponding generating function 
Q\(xi) and second-order cumulant «2 ( n= l) = af, in terms of which the generating function 
for x', the dual to q' = ^ Qi/(,o'iy/n), is 



1=1 



= n / ex p 



iqiX' 



(Tly/n 



(25) 



For such convolutions, cumulants of f(q') are related to cumulants of f(qi) by K^\n) = n«' 9) (l), 



so that 



0(1) 



n2" 



(26) 
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and hence depends on £ = n 1 ^ 2 through 



<S>(x') 



-x' 2 /2 



cxp 



3=3 



e x /z exp 



i=3 



(27) 



Expanding the exponential in powers of £ rather than of x' and again inverting term by term, 
we obtain the Gauss-Edgeworth series 



/(</) = foW) 



1 + £ 2 QH 4 + £ 4 (±$H 8 + ( e H 6 ) + ^ 6 (Kl^ + UC 6 H W + Cs^s) 



+£ 8 {^CiHiQ + \ciC&Hi4, + + C4C8-ffi2 + Cio-ffio) + • • • 



(28) 



again writing H r = H r (q') for short and with Q here understood as Cj(l)- Unlike the equivalent 
GGC expansion of Eq. (|16p . in which the order of the expansion was determined by the order of 
H m , a given term of order £ w in the GEW series is a linear combination of Hermite polynomials 
of different order. 

The relation between Gram-Charlier and Edgeworth ordering is summarised in Table 1, 
with terms listed in ascending order for w. A given term CriCr 2 - " " is characterised by the set of 
partition coefficients = 4, 6, 8, . . . which are constrained to the Gram-Charlier and Edgeworth 
orders by 



k 

k 



m . 



w . 



The re-ordering becomes important already for the second-lowest order w = 4. 



(29) 
(30) 



series 




Edgeworth 


Gram-Charlier 


term 


{r k } 


order w 


order m 


U 


{4} 


2 


4 


Ce 


{6} 


4 


6 


t 2 


{4,4} 


4 


8 


Cs 


{8} 


6 


8 


CeC 4 


{6,4} 


6 


10 


t z 
U 


{4,4,4} 


6 


12 


ClO 


{10} 


8 


10 


C 2 


{6,6} 


8 


12 




{8,4} 


8 


12 


CeCl 


{6,4,4} 


8 


14 


cl 


{4,4,4,4} 


8 


16 



Table 1: Re-ordering of terms between Gram-Charlier (GC) and Edgeworth (EW) series 
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4.2 Test using SNIG 

Edgeworth re-ordering of terms in the derivatives <S>v)(q f ) leads to expressions for the a>cumulants 
as ratios of power series^) in I. For the SNIG test case, these series simplify to 



(x) 
74 



i + h^ 2 + i# 



3072 U l 



+ 



_385_ f, 
98304 '4* 



+ 



1 + hd 2 



5 7 4 2 ^ 4 + gf^ 6 



175 



384 



98304 



Id 2 + + §§ 7 Jj 



35 „,4 
1152 



+ ••• 



1 + 5 /2 , 55^2/4 , ^5.^3/6 i 35 4^8 + 
4 ' ^96 '4* ^ 384 '4 C ^ 18432 '4* ' • 



(31) 
(32) 



with 74 = 7 4 9) (1). The algebraic simplicity of the above compared to the equivalent GGC 



relations pTp - (jl8j) and the GEW relation ([28]) is due to the fact that SNIG kurtoses obey 



(33) 



with constants {F 4 , F 6 , F 8 , . . .} = {1, 5, . . .} fully determined by the SNIG pdf. While rela- 



tion (|33p is of course fulfilled by all convolutions through (|26p . it is true for the SNIG case even 
without convolution. 

In Fig. 2, we show the percentage deviations of the Edgeworth-truncated approximations 
(I3ip - (|32p from their respective exact SNIG values as a function of the g-kurtosis . The 
improvement in accuracy over the GGC ordering is dramatic. Unlike the GGC, the GEW series 
also continues to improve as higher orders of w are included. 



4.3 n-divisibility and the Central Limit Theorem 

What structure or principle underlies the strong superiority of GEW over GGC ordering? 
Clearly, the expansion parameter n = £~ 2 must be playing a crucial role. Eq. (125p characterises 
the generating function ^(x') and therefore the pdf f(q') as "n-divisible" or more prosaically as 
precisely that convolution f(q') = 5[q' — J2i ( lil \\pn\ Ili f(li) which led to Eq. ([25]) . The GEW 
method might therefore be expected to work well if n-divisibility could be established for a given 
experimental data set. However, it is usually not possible to directly establish whether the f(q') 
of an experimental data set is n-divisible, and there is no obvious physics reason to believe that 
a momentum difference q' between two particles is the result of an underlying summation. 

The reason for the success of GEW relies not on physics assumptions, but on a better 
description of nongaussian systems. In statistics terms, any deviation from the gaussian is 
captured by the GEW series in the rate of approach of higher-order cumulants to zero with 
increasing n. Indeed, many proofs of the Central Limit Theorem rely on the fact that Kj(n) = 
nKj(l) immediately yields the fixed rate of convergence for any generalised kurtosis as shown 
in Eq. ([26D . 

The success of GEW ordering is therefore based on the fact that all contributing terms in a 
given order of I = rT 1 ! 2 have the same rate of convergence to the gaussian limit. Furthermore, 
due to the alternating sign of i?2r(0) = (— l) r (2r — 1)!!, the sum of contributions within a given 

2 Since k± and contain products of generating functions, terms of order higher than £ w are generated. 
Such terms must of course be omitted in a consistent 0(£ w ) calculation. 
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Figure 2: Gauss Edgeworth series: Percentage deviations of GEW approximations of found 
by Eq. (J3T]) vs. (22]) and of 7^ found by d32j) vs. (|24p . as a function of the g-kurtosis 7^ for 
various Edgeworth orders w. The dashed line at zero marks perfect agreement between the 
truncated expansion and the exact SNIG answer. 



0(£ w ) term tends to be substantially smaller than the individual contributions; for example the 
w = 4 term for the SNIG test case is made up of ±Cf#8(0) = 0.09l7f and C6^6#6(0) = -0.1047f , 
adding up to —0.01 374. 

It is not necessary to know the value of n to make use of GEW ordering: Using (|26p , we can 
merely re-absorb the l w in (|28p to write the GEW expansion in n-independent notation (now 
with Q = Cj( n ), the measured kurtosis) 



fW) = M) 



1 + (C4# 4 ) + (K1^8 + ( G H 6 ) + (±<f #12 + C4C6#10 + C 8 H 8 ) 



+ (^^4-^16 + 5C4C6-H14 + + C4C8-£fl2 + ClO-fflo) + • • • 



(34) 



but keeping in mind that each expression within a round bracket represents a certain rate of 
convergence and must be included or excluded as a whole. 

It is not even necessary to require n-divisibility as such: For the GEW ordering to be effective 
we require only that f(q) is reasonably close to a gaussian, where "reasonable" is typically 
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quantified by the errors shown in Fig. 2. The derivation also does not rely on a particular form 
of /i(<?i) other than requiring existence of its cumulants. 

5 Conclusions 

We have calculated expressions for x-cumulants of second and fourth order for the emission 
region in terms of measured g-cumulants. On using a nongaussian test function to quantify 
accuracy of expansions, we have shown that the textbook Gram-Charlier series is unsuitable at 
any level of approximation. By contrast, the Gauss-Edgeworth expansion, which orders terms 
based on the rate of approach to a gaussian, does give results which become increasingly accurate 
as more terms are added. The GEW series is robust in all respects; for example it does not 
require n-divisibility as such but only that f(q) be close enough to a gaussian in the sense that 
the measured g-cumulants entering truncated GEW expansion should have reasonable error 
bars. 

The present one-dimensional calculation clearly cannot be applied immediately to experi- 
mental data, but is meant to show that, even on the fundamental level of expansions, there 
are major questions which must be addressed first. In sorting out the fundamental issue of 
re-ordering, the present results represent an important step towards a consistent framework for 
shape description. 

Application to experimental data will require generalisation to three dimensions using the 
existing 3D machinery of Refs |15} I18j . Furthermore, sampling fluctuations of experimental 
cumulants will have to be taken into account. In this connection, we also note that the GEW 
ordering has the additional advantage of placing terms with higher powers of 7^ into lower 
orders of w, making it unnecessary to measure higher-order kurtoses. Based on slightly different 
arguments, Cramer [20] also concluded that GEW was superior to GGC ordering; this has also 
been verified GGC vs GEW comparisons of the nongaussian pdf itself [214 122] , The Gram- 
Charlier series can therefore be considered to be inferior to the Edgeworth equivalent in all 
respects. 

Note that it is not necessary to measure the correlation function f(q) at q = 0, despite the 
fact that x-cumulants rely formally on the generating function fjllH at zero. The q-cumulants 
themselves are functions of f(q) over the whole range of q, while the correlation strength pa- 
rameter A cancels in the normalisation ([3]). 

Our final comment pertains to the usual practice of obtaining information on the correla- 
tion function through fits of nongaussian parametrisations. Fits rely on an a priori choice of 
parametrisation, guided only by the minimisation of x 2 > an d suffer from increasing ambiguity 
in higher dimensions. By relying on direct measurement of coefficients, the present method and 
those of Refs [111 [T2] etc leave less room for arbitrary choices and put the uncertainty where it 
belongs: in the sampling fluctuations of measured experimental quantities. 

Acknowledgements: This work was supported in part by the National Research Founda- 
tion of South Africa. 
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